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NUMERICAL SIMULATION OL ION THRUSTER OPTICS 


Cody C. Farnell*, John D. Williams, and Paul J. Wilbur 
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ABSTRACT 

A three-dimensional simulation code (ffx) designed to analyze ion thruster optics is described. It is 
an extension of an earlier code and includes special features like the ability to model a wide range of grid 
geometries including cusp details and mis-aligned aperture pairs. The motivation for advancing the code 
was the need to accurately predict grid lifetimes and to study the physical evolution of the surfaces that are 
being eroded and surfaces where net material buildup occurs. The techniques used to model the sputter 
erosion processes and track the sputtered material are described. The ability of the ffx code to track material 
eroded from the grids that flows downstream into the ion thruster plume would be valuable to interpret data 
collected using advanced laser diagnostic systems. Ground-based life testing of ion thruster optics, essential 
to the understanding of the processes of grid erosion, can be both time-consuming and costly. Often these 
tests only provide limited data that is collected at the end of the test when it is too late to make any necessary 
changes to flight hardware. Even worse, data on only one test article are typically obtained, which are often 
viewed suspiciously due to test facility effects and the obvious lack of statistical significance. This situation 
is exacerbated for life tests of future ion thruster optics that are being designed for very ambitious 
applications requiring much larger propellant throughput. The goal of the ffx development effort is to 
provide a simulation tool to the ion propulsion community that can be used in combination with short-term 
tests of ion thrusters that are designed to alleviate the need to perform full-blown life tests. In this paper 
results of simulations are presented that describe wear profiles for several standard and non-standard 
aperture geometries, such as those grid sets with square- or slotted-hole layout patterns. The goal of this 
paper will be to introduce the methods employed in the ffx code and to briefly demonstrate their use. 

THE FFX CODE 

Many simulation codes have been developed to study various aspects of ion thruster optics. One 
such code is the igx code, developed by Nakayama and Wilbur for the high-speed, three-dimensional 
analysis of grid sets with axially aligned, hexagonal aperture layouts. 1 This code has been shown to agree 
well with experimental data in high specific impulse applications. 2 The ffx code analyzes a three- 
dimensional, rectangular region with symmetry conditions applied on all sides. A uniformly spaced 
Cartesian mesh, efficient for cell indexing, 3 is applied to the volume with each direction having a different 
mesh spacing. For the simulation of a typical aperture pair in a hexagonal aperture layout, a mesh of 
approximately 30 by 50 by 300 cells is applied to model two quarter-sized apertures. A cross-section of the 
typical geometry and applied potentials for a two-grid system is shown in Figure 1. The key simulation 
tasks in the ffx program are outlined in the flowchart in Figure 2. The following discussion will elaborate 
upon various aspects of these tasks. 
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Figure 1. Typical two-grid thruster geometry. 



Figure 2. Flowchart of the ffx simulation code. 


For the case of an electrostatic thruster, where magnetic fields can be neglected in comparison with 
the electric fields in the region near the accelerating grids, mesh potential values, denoted by ()>, are described 
by the Poisson equation 


'V = -— 


or 


d 2 </> d 2 </> d 2 (b 

+ — V + - 


dx 1 


dz 2 


p_ 


( 1 ) 


-o dy kja o 0 

In this equation, p is space charge and So is the permittivity of free space. In the ffx code, both singly 
and doubly charged ions as well as electrons will contribute to the overall space charge. Thus the right-hand 
side of this equation will be 


P _ P+ + P++ + Pe , . 

VW 

where p+ and p++ are the space charge contributions from singly and doubly charged ions respectively, and 
p e is the space charge contribution from electrons. 


Poisson’s equation is solved for mesh point potentials at the start of each beamlet loop. Initially, no 
ion or electron space charge is applied upon the grid mesh. Consequently, the potential solution routine of 
the code first solves the Laplace equation, where the space charge is zero. 


NASA/CR— 2003-212305 


2 


Electric fields in the region are related to potential by 
E = -V 0 or E = —| 


r d(j) ? d(j) - 3^ 
— z +— i + —k 
dx dy dz j 


( 3 ) 


The second order partial derivatives in the Poisson equation and the first order partial derivatives in the 
electric field equation are approximated by finite difference equations of the central difference type. 4 For 
instance, in the x direction: 
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Following other simulation codes, 3 ' 5 equations relating the electron space charge to ion space charge 
at each mesh point are solved within the potential solution routine of the code depending on the average 
upstream and downstream values of space charge and the current value of each mesh potential. In the region 
upstream of the grids within the discharge chamber plasma: 
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In the region downstream of the grids within the beam plasma: 
Pe =-(p + oo +P ++ oo)exp 
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In these equations, <|)o is the value of the upstream plasma potential, (j)* is the value of the 
downstream plasma potential, <|) is the current value of the mesh potential, T e0 is the upstream electron 
temperature, and T eoo is the downstream electron temperature. 


The Gauss-Seidel method is used to solve Poisson’s equation. In this method, Poisson’s equation is 
rearranged to obtain (jy^ once the finite difference approximations for the partial derivatives have been 
substituted. Iterations are then performed through all mesh points that are not at a known potential until all 
of the mesh point potentials are changing by a value less than a preset limit. Relaxation is used to assist 
convergence of the solution by intuitively predicting each mesh point’s next potential value based on the 
current and last potential values. Relaxation is implemented using 

tnew =«^ + ( 1 -«V 0 H Where 0 <a <2 (7) 

In this equation, a is the called the relaxation parameter. Without the electron density equations embedded 
within the potential solution, as in the first iteration, using a value for a between 1.0 and 2.0 will increase 
the rate of convergence, called overrelaxation. With the electron density equations added however, the 
electron density will change along with mesh potentials during the potential solution, making the solution of 
the equation non-linear. In this case, a value of a that is between 0.0 and 1.0 is used for convergence, called 
underrelaxation, because it helps to damp out oscillations that the electron equations create. 


Singly and doubly charged ions are represented by particles that are injected into the upstream 
boundary and terminate upon one of the grids or the downstream boundary. During the movement of these 
particles, they have the mass and charge of a single ion, whether they are singly or doubly charged. Each of 
these particles carries with it a much larger charge than a single ion in order to greatly reduce the number of 
particles required to simulate the ion beamlet. For a typical simulation, about 15,000 particle trajectories are 
tracked within each beamlet loop. Particles are injected into the upstream boundary of the analysis volume 
with axial velocities equal to the Bohm velocity: 
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where q is the particle’s charge, kT e o is the upstream electron temperature, and m is the mass of the particle. 


Forces on charged particles are described by the Lorentz equation: 

F = q(vxB + E) (9) 

In this equation, q is particle charge, v is particle velocity, B is the magnetic field, and E is the electric field. 
In the electrostatic case, again neglecting magnetic field effects, the vector equations governing the motion 
of particles are given in equations 10a and 10b. The resulting equations of motion in the x direction, for 
example, are given in equations 1 la through lie. 


F = qE = m — 
dt 



(10a) 

(10b) 
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The time step At is constantly adjusted during particle motion to ensure that a particle’s movement is 
kept under a specified fraction of the smallest mesh spacing, h, in accordance with the Courant condition, i.e. 


Av • At < 



( 12 ) 


Making sure particles do not travel beyond a certain distance in any given time step is required to ensure that 
details of the system are not lost. 


Throughout particle movement, a particle’s total energy remains constant and is given by the sum of 
its kinetic and potential energies: 


E = —mv 2 + qd) 

2 


(13) 


Using the preceding equations of motion, however, the total energy of a particle is not guaranteed to remain 
constant. Thus at the conclusion of each time step, its total energy is forced to be constant by adjusting the 
length of its velocity vector according to the particle’s total initial energy and local value of potential. 
Specifically, the particle’s three velocity components are multiplied by the ratio of the required velocity 
magnitude for energy conservation to its current velocity magnitude in order to maintain the same direction 
of motion while keeping the total energy of the particle constant. 


As particle trajectories are followed through cells, the charge that each particle is carrying is applied 
to the eight surrounding mesh points by the method of volume weighting. In this method, the charge applied 
to each mesh point is proportional to the fraction of the total cell volume that the volume closest to the 
opposite comer occupies. For instance, the charge applied to Point 1 in Figure 3 is proportional to the 
fraction of the total volume that the volume near Point 8 occupies: 


C ^§ Point 1 


C ^S Particle 


(x 2 -x)-(v 2 -y)-(z 2 - z ) 
Av • A y ■ Az 


(14) 


During particle movement, the electric fields in each direction at the particle’s actual location are 
interpolated from the known electric fields at the eight comer mesh points of the cell using the same 
weighting scheme that is used to apply charges onto the mesh points in a reverse manner. The electric field 
at the particle in the x direction, for example, is the sum of the eight x direction electric fields at the comer 
points multiplied by the eight respective volume fractions. 


The neutral density in each cell, n n , is calculated by assuming a radially uniform, free molecular 
flow of neutral particles. Clausing factors, which are dependent on each grid’s thickness and aperture area, 
are used to adjust the neutral density through the grids. An upstream neutral density can be specified, or the 
beamlet current and propellant utilization efficiency can be used to determine the flow rate through the 
region. 
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Once singly and doubly charged ion densities and the neutral density in each cell have been 
calculated, equations describing the volumetric charge exchange ion production rates in each cell are used to 
introduce charge exchange ions. Charge exchange ions are created when electrons transfer to fast moving 
ions from relatively slow moving neutral atoms, creating slowly moving ions and fast moving neutral atoms. 
Charge exchange ion production rates are described using equations of the form 


dt 

Depending on the reaction in question, n, can be the singly or doubly charged ion density, Vj the velocity that 
a singly or doubly charged ion would be moving through the cell, and g(v ;) the cross section of the specific 
reaction. 

Erosion of the grids can be due to both ions from the upstream plasma as well as charge exchange 
ions impacting the grids. The number of sputtered grid atoms is dependent on the energy and angle of 
incidence of the impacting particle. The angle of incidence is found using the local surface normal vector of 
the impacted cell. Because erosion will constantly change the grid surface, surface normal vectors are found 
at any time by calculating the location of the regional center of mass of the impacted cell relative to its 
geometrical center. The surface normal of the impacted cell is then defined as the line that extends from the 
center of mass in the region through the center of the impacted cell, as shown in Figure 4. 




(v,) 


(15) 



Figure 4. Surface normal vectors and angle of incidence. 

Grid atoms are sputtered with a cosine distribution based on the surface normal vector of the impacted cell 
and are subsequently followed from the particle impact point until they deposit onto one of the grids or exit 
through the upstream or downstream boundaries. 


NASA/CR— 2003-212305 


5 


Charge exchange ions originating from several centimeters downstream of the grids can flow back 
to the accel grid and cause erosion. To keep an analysis volume of reasonable size and a simulation of 
reasonable speed while still simulating charge exchange ions from far downstream, charge exchange ion 
production is adjusted near the downstream boundary to bring the current reaching the accel grid to the level 
measured during experimental tests conducted in vacuum facilities with very high pumping speeds. 

RESULTS 

To illustrate the use of the ffx code, simulations were run on four SUNSTAR (Scaled Up NSTAR) 
grid sets 6 that have various aperture patterns. The geometries of these grids and the simulated operating 
conditions are shown in Figure 5. The parameter chosen to compare the grid sets is normalized perveance 
per unit grid area: 



In these equations, J b is the beamlet current, A TG the aperture and web area, l e the effective grid spacing, and 
V T the total accelerating voltage. A comparison of the beamlet shapes at the perveance and crossover 
impingement limits is shown in Figure 6. At the perveance limit, where extraction currents are relatively 
high, ions impinge upon the accel grid on the same side of the aperture from which they started. At the 
crossover limit, where extraction currents are relatively low, ions can crossover the aperture centerline and 
impinge upon the opposite side of the accel grid. At the center of a thruster, aperture pairs would likely be 
operating closer to the perveance limit where ion densities are high. Conversely, at the edges of the thruster, 
aperture pairs would be operating closer to the crossover limit where ion densities are lower. 

Figure 7 compares the impingement limits found experimentally at CSU 6 with those found using the 
ffx code for the four SUNSTAR aperture patterns. The ffx code generally predicted slightly greater 
perveance and crossover limits than those found experimentally. Overall, good agreement was seen between 
the relative limits among different aperture patterns. 

Figure 8 compares the sheath and neutralization surface shapes found near the perveance and 
crossover limits for each of the aperture patterns. The sheath and neutralization surfaces shown here are 
surfaces of constant potential: the sheath at the upstream plasma potential and the neutralization surface at 
the downstream plasma potential. 

Erosion tests were run on the four SUNSTAR grid sets at beamlet current values near the perveance 
limit for each grid set. In each case, the current reaching the accel grid was adjusted to 0.3% of the beamlet 
current. Figure 9 compares the erosion patterns produced on each of the grid sets at two values of propellant 
throughput, the second being twice as much as the first. Note that each of the grid sets would be run for 
different amounts of time to reach the same level of propellant throughput. Because the grids are operating 
near the perveance limit, high-energy charge exchange ions produced in the region between the grids quickly 
erode away the accel grid barrel from the upstream side of the accel grid in all of the aperture layouts. 
Following this period of immediate barrel erosion, charge exchange ions from the downstream region create 
erosion patterns in the downstream side of the accel grid with less accel grid barrel erosion. The familiar pit 
and groove pattern of erosion can be seen for the hexagonal aperture pattern. 

CONCLUSIONS 

A three-dimensional optics simulation code was developed that can be used to simulate ion optical 
behavior in complex grid geometries with reasonable efficiency. Results from the code were found to be 
consistent with experimental data in predicting impingement limits. Erosion patterns predicted by the code 
agree well with the known erosion patterns of grids with hexagonal aperture layouts, and the code yields 
erosion patterns for grids with square and slotted aperture layouts that are judged to be reasonable. 
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Full Slot Pattern 



Geometry / 
Operating Conditions 

I = 5000 s 
V p = 2.91 kV 
V s = 2.88 kV 
V a = -360 V 
V d = 0V 
d s = 3.1 mm 
d a = 1.86 mm 
t s = 0.62 mm 
t a = 0.83 mm 
£ cc = 3.6 mm 

Vif =4 - 8mm 

= 12.0 mm 


l j = ffx analysis region 

I I = region displayed in figures 

Figure 5. SUNSTAR grid aperture patterns. 



Beamlet shape at the 
perveance limit 


Beamlet shape at the 
crossover limit 


Figure 6. Perveance and crossover impingement limits. 
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IMPINGEMENT TO BEAMLET CURRENT RATIO [J imp /J b ] (%) IMPINGEMENT TO BEAMLET CURRENT RATIO [J imp /J b ] (%) 


Experimental 



NORMALIZED PERVEANCE PER UNIT GRID AREA [P TG ] 

ffx Code 



NORMALIZED PERVEANCE PER UNIT GRID AREA [P TG ] 


Note: Charge exchange impingement current was not included in the ffx data. 
Figure 7. Comparison of experimental and ffx code predictions of crossover and perveance limits. 
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Crossover Limit 


Perveance Limit 


Hexagonal Pattern 


P TG = 0.13 
J b = 0.25 mA 


Square Pattern 


P TG = 0.13 
J b = 0.29 mA 



P TG = 0.31 
J b = 0.59 mA 


P TG = 0.26 
J b = 0.57 mA 



Sheath 

Screen Grid 

Accel Grid 

Neutralization 

Surface 





Figure 8. Sheath and neutralization surface shapes for the four SUNSTAR grid sets. 
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Propellant Usage 1 
(-5000 hours of operation 
for hexagonal grid) 


Propellant Usage 2 
(-10000 hours of operation 
for hexagonal grid) 



Upstream 

Cells 



Half Slot Pattern 

? TG = 0.32 
J b = 1.64 mA 




I 

Downstream 

Cells 



Note: Colored cells are cells that have been eroded away. 
Figure 9. Accel grid cell erosion for the four SUNSTAR grid sets. 
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